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Abstract 


The number of neuroimaging data sets publicly available is growing at fast rate. The increase in 
availability and resolution of neuroimaging data requires modern approaches to signal processing for 
data anal ysis and results validation. We introduce the application of sparse multiway decomposition 
methods iCMiafa and Cichockl, hoii) to linearized neuroimaging models. We show that decomposed 
models are more compact but as accurate as full models and can be successfully used for fast data 
analysis. We focus as example on a recent model for the evaluation of white matter connectomes 
jPestilli et aihoid) . We show that the multiway decomposed model achieves accuracy comparable to 
the full model, while requiring only a small fraction of the memory and compute time. The approach 
has implications for a majority of neuroimaging methods using linear approximations to measured 
signals. 

Keywords. Brain anatomy, Big data. Brain mapping. Brain networks, Connectomics, High-performance 
computing. Individualized medicine, Myeline, Personalized medicine. Precision medicine. Tensor de¬ 
composition. 


1. Introduction 

NEuroscience is transforming. Investiga¬ 
tors have recently started sharing brain 
data collec te d on la rge human populations 
i Scott etal. I2OIII: IVan Essen and Ugurbil . 


2OI2I: Van Essen et at 2OI3I) . The moderri 

era of data sharing has the potential to pro¬ 
mote scientific replicability of results and 
advance our understanding of brain func¬ 
tion. Shared data can be accessed by differ¬ 
ent research groups and analytic tools can 
be tested and applied to the very same data 
in atte mpt to extend and replicate sc ientific 
results I Pestillil boisl: l^o et al , 2014 ). This 
fundamental shift of modern neuroscience 
will allow separating the wheat from the 
chaff. This paradigm shift is changing the 
landscape of brain imaging and motivating 


investigators to implement new met hods for 
data storage, analysis and sharing (jPestilli , 
201.51) . 


Hereafter, we show that modern 
neuroimagmg can benefit from re- 
cent developments m sig n al processing 
llCaiafa and Cicho^ . 120121: ICichocki et al , 
1201.51) to address the new challenges raised 
by the most recent growth in data size, reso¬ 
lution and neuroscience models complexity. 
Standard signal processing methods in neu- 
roimaging are based on classical Linear Al¬ 
gebra and use mathematical objects such as 
one dimensional vectors (ID) and two dimen¬ 
sional matrices (2D). Multiway approaches 
exploit the properties of arrays of higher or¬ 
der; arrays with three or more dimensions. 
Multiway arrays, also referred to as hyper¬ 
matrices or tensors, are the generalization 
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■ Thalamic radiation ilFOF DCingulum cingulate ■ SLF _ 

□ Forcepts {M and m) ilLF iCingulum Hippocam ■ Uncinate 10 50 90 

■ Corticospinal tract ■ Arcuate Probabilistic 


Figure 1: Tractography comparison and evaluation. A. Ten major human white matter tracts gmerated using 
probabilistic tractography. B. Ten major human white matter tracts generated using deterministic tractography. C. 
Comparison ofr.m.s. errors in predicting the diffusion signal of probabilistic and deterministic connectomes. 


of vectors and matrices to higher number 
of dimensions (ND). Multiway arrays have 
two primary advantages when applied to 
neuroimaging data compared to classical lin¬ 
ear algebra approaches. First, neuroimag¬ 
ing data are multidimensional in nature and 
their structure is preserved by multiway ar¬ 
rays, which in turns allows for straightfor¬ 
ward data addressing. Second, they pro¬ 
vide convenient and compact representations 
of neuroimaging models that allows sub- 
stantial reduction in memo ry consumption 
I Caiafa and Cichocki , 2013t Cichocki et al . 
2015 : Morupl, 2011 ). Below we briefly in¬ 


troduce basic concepts for multyway array 
decomposition and show one example ap¬ 
plication of the methods to mo deling diffu¬ 
sion imaging and tractography dPestilli et al 


2014). 


dMRl data in combination with fiber 
tracking allows measuring the anatomy 
and tissue properties o f the human white- 
matter in liv i ng brains dYeatman et al , l2012t 
Yendiki et all. I2OIII: Izhang et all, 2010 . By 


measuring living brains this technology 
allows correlating the properties of the 
white matter tissue and structure with hu¬ 
man behavior, cognition as well as de- 
velop ment and aging in health and dis¬ 
ease d Sa manez-Larkin an d Kn uts^. 20151: 
Thomason and Thompson ! l2011 ). dMRl and 


tractography generate estimates of white 
matter tracts and connections, such estimates 
need me thods for routine evaluation and val¬ 
idation I Catani et al , 20121: lones et al . 20131: 
Sporns et aj 2OO5I) . For example, it is estab¬ 


lished that different tractography methods 
can generate different estimates of the shape 


of the white matter anatomy (Fig. [T]4 and B). 


Recently we developed a m ethod called 
LIFE, Linear Fascicle Evaluation dPestilli et al , 
201 4 . LiPE is an approach to statistical val¬ 
idation of in vivo connectomes that can be 
applied to individual brains. To validate 
connectomes in each brain, LiFE requires as 
input a set of candidate white-matter fasci¬ 
cles genera ted using available tractography 


algorithms llCook et al 

0 

0 

liane et ai 

2006; 

Tournier et ai 12012). 

The candidate fasci- 


cles are used to generate a prediction of the 
anisotropic diffusion signal measured within 
the white matter volume using dMRl. The er¬ 
ror (root mean squared error, r.m.s.) in pre¬ 
dicting the diffusion signal is computed us¬ 
ing cross-validation and used to establish the 
accuracy of a tractography solution. Por ex¬ 
ample, Pig. [T]L shows a comparison of the 
cross-validated r.m.s. error of two tractogra¬ 
phy methods in the same brain (the models 
shown in Fig. [l]4 and[Tp). The result demon¬ 
strate that in a majority of the white matter 
the model in Fig. [T]4 has lower r.m.s. error 
in predicting the diffusion signal. 


LiFE predicts the demeaned dMRl data 
by finding the linear combination of the dif¬ 
fusion prediction contributed by all individ¬ 
ual fascicles within a connectome. 
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Fascicle predition (diffusion tensor) 
Fascicle orientation 
a, /3 : Spherical coordinates 


Voxel V 



Demeaned signal 
at voxel V 
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Contributions of 
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C Demeaned 
diffusion signal 

[NeN^ X 1) 


Tableau Matrix M 

(iV„]V, X N.) 

h h 



■ ■■ Fascicle signal contributions. 

■ Demeaned signal (measurment) 
□ Empty entries (zero values). 


D 50 


CO 


10 


• • 


• Subject 1 

• Subject 2 

• Subject 3 

• Subject 4 

• Subject 5 


STN96 STN150 


Figure 2: The linear fascicle evaluation method (LiFE). 
A. Idealized diffusion prediction model for the fas¬ 
cicles in a voxel. The fiber-predicted diffusion sig¬ 
nal in a voxel is mod eled by the Steiikal-Tanne r 
equation (orange surface; Basser and PierpaoJl Il201l\) : 
\Steiskal and Tannen ill 9651) ). evaluated in the direction 
of each fascicle in the voxel (red arrow). B. The dif¬ 
fusion prediction in a voxel with multiple fascicles is 
computed as a linear combination of the prediction of 
each individual fascicle (Eq. O. The demeaned diffu¬ 
sion signal (brown) is predicted as the weighted sum 
of the predictions from all fascicles (green and red) in 
the voxel. C. Global linear model. The demeaned mea¬ 
sured diffusion signal across all white-matter voxels is 
organized in a long ID vector (y). The signal predic¬ 
tion from all fascicle across the full white-matter vol¬ 
ume is organized into a large 2D block-sparse matrix 
M. Weights are assigned to each fascicle using an opti¬ 
mization procedure, fascicle weights are organized into 
a ID vector (w). Notations, Ng represents the number 
of diffusion directions, Nj, the number of voxels, and 
Nf the number of fascicles). D. LiFE storage measure¬ 
ments. The size of the LiFE Model (2D matrix M in 
Panel C) measured in GB in several brains and data 
sets. Matrices built using double floating-point preci¬ 
sion. 


The fundamental problem solved by LiFE 
is to allow computing an error measure 
that can be used for a variety of tasks in 
evaluating individual cormectomes proper- 


ties generated with virtually any tractogra- 

ohy method dComez et al. 

2015 

: Pestilli et al. 

2014 

; iTakemura et al. 20151 ' 

ieatman et al. 

2014 

). One consequence of formulating the 


LiFE model as a system of linear equations 
is the size of the 2D matrix representing the 
model (see Fig. |2]d and Eq. |8j. The compu¬ 
tational demands of the LiFE method repre¬ 
sented as 2D matrix are substantial. Given an 
approximate human white-matter volume of 
500m/, the spatial and directional resolution 
of modern dMRI measurements and a rea¬ 
sonable number of fascicles in a cormectome, 
the size of the LiFE model can be as large as 
about 50GB (see Fig. |2)D). 

Hereafter, we introduce a multiway de¬ 
composition method that reduces the size of 
the LiFE model by over 97%. The new multi¬ 
way approach maintains memory consump¬ 
tion under 1GB per brain. The approach 
can be applied to any linearized model for 
global tractography, tractography evaluation 
and microstructure estimation. Below, we (1) 
describe the mathematics of the sparse mul¬ 
tiway decomposition, (2) show that memory 
consumption is quasi-constant as function of 
connectome size when using the decomposi¬ 
tion and (3) replicate major results of the orig¬ 
inal LiFE work using the new appr oach and 
several datasets dPestilli et all. I20l3) . We pro¬ 
vide open source software implementing the 
LiFE model at github.com/brain-life/life 
and f rancopestilli. github. io/lif e Scripts 
and data used to generate the results in the 
present articles can be found at these reposi¬ 
tories. 


II. Materials and Methods 

I. Diffusion-weighted MRI acquisi¬ 
tion 

Diffusion-weighted Magnetic Resonance 
Imaging data (dMRI) were collected in 
five males subjects (age 37-39) at the 
Stanford Center for Cognitive and Neu- 
robiological Imaging using a 3T Gen- 
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eral Electric Discovery 750 (General Elec¬ 
tric Healthcare) equipped with a 32- 
channel head coil (No va Medical). This 
dataset is publi shed (jPestilli et al , 20141: 


Rokem et al 1201511 and publicly available at 
http://purl.stanford.edu/rt034xr8593 and 
http: //purl. stanford.edu/ng782rw8378 , Data 
collection procedures were approved by 
the Stanford University Institutional Review 
Board. Written consent was collected from 
each participant. 

Stanford 96 diffusion directions data set 
(STN96): two diffusion-weighted with whole- 
brain volume coverage were acquired in five 
individuals within a single scan session us¬ 
ing a dual-spin echo diffusion-weighted se¬ 
quence. Water-proton diffusion was mea¬ 
sured using 96 directions chose n using the 
elect rostatic repulsion algorithm ( tones et~^. 


19991) . Diffusion-weighting gradient strength 
was set to 2, OOOs/mm^ (TE = 96.8ms). Data 
were acquired at l.Snim^ isotropic spatial res¬ 
olution. Individual data sets were acquired 
twice and averaged in k-space {NEX = 2). 
Ten non-diffusion-weighted (b = 0) images 
were acquired at the beginning of each scan. 


Stanford 150 diffusion directions data set 
(STN150): for one subject two dMRI data sets 
were also acquired in a single session using 
150 directions, 2mm^ isotropic spatial reso¬ 
lution and b values of 2,000s/nirff {TE = 
83.1,93.6, and 106.9ms). 


MRI images for STN96 and STN150 were 
corrected for spatial distortions due to BO 
field inhomogeneity. To do so, BO magnetic 
field maps were collected with the same slice 
prescription as the dMRI data using a 16- 
shot, gradient-echo spiral-trajectory pulse se¬ 
quence. Two volumes were acquired before 
and after the dMRI scan {TE = 9.091ms and 
11.363ms respectively), the phase difference 
between the volumes was used as an esti¬ 
mate of the magnetic field. 

Subjects' motion was corrected us- 
ing a rigi d -body alignment algorithm 
jEriston et all l2004h . Diffusion gradients 
were adjusted to account for the rotation 
applied to the measurements during mo¬ 
tion correction. The dual-spin echo se¬ 
quence we used does not require eddy cur¬ 
rent correction because it has a relatively 


long delay between the RE excitation pulse 
and image acquisition. This allows for suf¬ 
ficient time for the eddy currents to de¬ 
phase. Processing software is available at 
https://github.com/vistalab/vistasoft 


II. Anatomical MRI acquisition and 
tissue segmentation 


The white/gray matter border was defined 
on the average of two 0.7-mm3 Tl-weighted 
ESPGR images acquired in the same scan 
session. Tissue segmentation was performed 
using an auto mated procedure (EreeSurfer 
Eischl 1 I 2 OI 2 1 I and refined manually 


(http: //www. itksnap.org/pmwiki/pmwiki.php). 


III. Whole-brain coimectomes gen¬ 
eration and visualization 


Eiber tracking was performed using the MR- 
trix 0.2 toolbox ( Tournier et al . 2012 ). White- 
matter tissue was identified from the cor¬ 
tical segmentation performed on the Tl- 
weighted images and resampled at the res¬ 
olution of the dMRI data. Only white- 
matter voxels were used to seed fiber 
tracking. We used two tracking meth¬ 
ods: (i ) tensor-bas e d deterministic t r actog- 
raphv llBasser et a\ . 2000t Ihazar et a 1 l200i. 


Tournier et al . 2012h and (ii) CSD-based 


probabilistic tracking ([Behrens et al. 


Parker et al , 2003 ; Tournier et all. 2012h with a 


maximum harmonic order of 10 {E^ax = 10, 
step size: 0.2mm; minimum radius of curva¬ 
ture, 1mm; maximum length, 200mm; mini¬ 
mum length, 10mm; fibers orientation distri¬ 
bution function {four) amplitude cutoff, was 
set to 0.1). 


We created candidate whole-brain 
connectomes with 500,000 fascicles in 
each individual brain (five), data set 
(two) and tractography method (two). 
Analysis were performed independently 
for each brain. Pigures of tracts 
and brain images were generated us¬ 
ing the Matlab Brain Anatomy toolbox: 
https://github.com/francopestilli/mba 
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IV. Brief introduction to modeling 
magnetic resonance diffusion signals 
from living human brain tissue 

The brain tissue comprises different cell 
types (e.g., neurons, astrocytes and oligo¬ 
dendrocytes). dMRI measures signals that 
depends on the combination of all cellular 
componenfs within the brain tissue. The 
dMRI measurements are generally modeled 
as the linear combination of fwo primary 
componenfs. One componenf describes fhe 
directional diffusion signal and is presum¬ 
ably relafed primarily fo the direction of the 
neuronal axons wrapped in myelin sheaths 
(white matter). This signal is often re¬ 
ferred to as anisotropic diffusion. The other 
component describes isotropic diffusion (non- 
directional) and is presumably related to the 
combination of signals originating from the 
rest of fhe cellular bodies within the brain 
tissue. Below we describe the simplest, fun¬ 
damental equations used to model the mea¬ 
sured dMRI signal in relation to the brain tis¬ 
sue components. 

dMRI measures the diffusion signal with 
and without diffusion sensitization. The 
measured signal results from fhe combina¬ 
tion of diffusion gradienf strength, duration 
and depends on the interval between two ap¬ 
plied gradient pulses. Below we denote the 
strength of sensitization with b and the vec¬ 
tors of measured diffusion directions (direc¬ 
tions along which the diffusion-sensitization 
gradients are applied) using the unit-norm 
vector 6 e E,^. 

For a given sensitization strength b and 
diffusion directions G, the diffusion signal 
measured at each location within a brain 
(voxel v) can b e computed using the fol¬ 
lowin g equation I Behrens et all. l2003bl: iFrank . 

2n02ll: 


S{9,v) fa woSo(w)e + X] be^Q/B^ 

fev 

( 1 ) 

where / is the index of fhe candidate white- 
matter fascicles within the voxel, So(u) is the 
non diffusion-weighted signal in voxel v and 
Aq is the isotropic apparent diffusion (diffu¬ 
sion in all directions). The value G^Q^G > 0 
gives us the apparent diffusion at direction G 


generated by fascicle f.Qf& is a sym- 
met ric and posit i ye-de finite matrix called ten¬ 
sor iBasser et al . 1994h . The tensor allows a 
compact representation of the diffusion sig¬ 
nal measured with dMRI. For example, Qf 
in the exponent of Eq. [T]can be replaced with 
the following simple tensor model: 



" Sa 0 O' 


’ Ui 

[ui U2 U3] 

0 Sri 0 


U2 


_ 0 0 Srj . 


. “3 . 


( 2 ) 

where Un G E^^^ are the unit-norm orthog¬ 
onal vectors that correspond to the semi¬ 
axes of the diffusion tensor ellipsoid, and Sa, 
Sr 2 define fhe axial and radial diffusiv- 
ity of the tensor, respectively. In the sim¬ 
plest version of the fascicle model, Sa = 1 
and Srj = Sr 2 = 0, which means thaf diffu- 
sion is restricfed to the main axis dire ction 
llBehrens et al , 2003b : Pestilli et al , 2014 ). 


V. Brief introduction to multiway ar¬ 
rays 

Multiway arrays generalize vectors (ID ar¬ 
ray) and matrices (2D array) to arrays of 
higher dimensions, three or more. Such ar¬ 
rays can be used to perform mulfidimen- 
sional facfor analysis and decomposition and 
are of interest to many scientific disciplines 
( Cichocki et allioiM Kolda and Badeii 120091) . 


Below we introduce a few basic concepts and 
notation helpful in discussing multiway ar¬ 
rays and their decomposition. We refer the 
reader to Table [T] for a summary of basic no- 
fations and definitions. 

Multiway arrays. Vectors (ID arrays) and 
matrices (2D arrays) are denoted below us¬ 
ing boldface lower- and upper-case letters, re¬ 
spectively. For example x G E^ and X 6 E^^l 
represent a vector and a matrix, respectively. 
A multiway array, also called N-th order ar¬ 
ray, is the generalization of matrices to more 
than two dimensions and is denoted by an 
rmderlined boldface capital letter, e.g. X G 
E^x/xK jg 2 _j.j order array of real num¬ 
bers. Elements {i,j,k) of a multiway array 
are referred to as Xjjj^. 

Array slices. Array slices are used to ad¬ 
dress a multiway array along a single dimen¬ 
sion (a cut through a single dimension of the 
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array). Slices are obtained by fixing the in¬ 
dex of one dimension of the multiway array 
while letting the other indices vary. For ex¬ 
ample, in a 3-rd order array X £ 
we can identify horizontal (f), lateral ()) and 
frontal (k) slices by holding fixed the corre¬ 
sponding index of each array dimension (see 

Fig.|3K). 

Unfolding arrays and mode-n vectors. Mul¬ 
tiway arrays can be conveniently addressed 
in any dimension by means of mode-« vec¬ 
tors. Given a multiway array X £ its 

mode-n vectors are obtained by holding all 
indices fixed except one, thus corresponding 
to columns (n = 1), rows (n = 2) and so on 
(see Fig. |^). For example, a 3D multiway 
array can be converted into a matrix by re¬ 
arranging its entries (unfolding). The mode- 
n unfolded matrix, denoted by X^^j £ 
where In = Ylm^n In and whose entry at row 
i„ and coluirm — 1 )J 2 ■ ■ • + 

-h (fw-i - 1)^N + In is equal to 

For example, mode-2 unfolding builds the 
matrix X( 2 ) where its columns are the mode- 
2 vectors of the multiway array and the rows 
are vectorized versions of the lateral slices, 
i.e. spanning dimensions with indices i and 
k (see Fig. |^). 

Multiway array by matrix product. Multi¬ 
way arrays can be multiplied by matrices (2D 
arrays) only if the matrix and multiway ar¬ 
ray sizes match in the mode specified for the 
multiplication. This is a generalization of ma¬ 
trix multiplication. Given a multiway array 
X £ ]Rhxh -x^N and a matrix A £ the 

mode-n product 

Y = X X„ A £ Rhx-x7„_ix/xl„+i-x7,v 
is defined by: 

yil-in-ljin+v-iN ^ (4) 

!„ = 1 

with 4 = 1/2,..., 4 ^ w) and j = 1,2,...,/. 

It is noted that this operation involves the 
products of matrix A by each one of the 
mode-n vectors of X- lo Fig. |^, the 3-rd 
order array by matrix product in mode-2 is 
illustrated, i.e. Y = X X 2 A. The size of the 
resulting array Y is that of X in mode-1 and 
-3 while the size of Y in mode-2 is /, i.e. it is 


equal to the number of rows in A. 



■ Horizontal slice 

■ Lateral slice 
IB Frontal slice 



□ Mode-2 vector 

□ Mode-3 vector 


L(2) 


Mode-2 

vector 



Figure 3: Multiway array decomposition basic nota¬ 
tion and operations. A. Frontal, lateral and horizontal 
slices of an example multiway array. B. Examples of 
mode-n vectors. C. Illustration of the mode-2 unfold¬ 
ing matrix X( 2 )' D- Array-by-matrix product (example 
product in mode-2). See Table\^or additional informa¬ 
tion about notation and mathematical definitions. 


III. Results 


I. The Linear Fascicle Evaluation 
model 


The Linear F a scicle Evaluation method (LiFE; 
Pestilli et all ( 2014 )1 allows evaluating the 
accuracy of connectomes generated using 
dMRl and any fiber tracking algorithm. The 
method evaluates the tractography solution 
by estimating the contribution to predicting 
the measured diffusion signal from all fasci¬ 
cles contained in a connectome (see Eig. |2]4 
and B for examples of fascicles). The eval¬ 
uation method focuses on the fascicles, for 
this reason it predicts the anisotropic diffu¬ 
sion signal; the demeaned diffusion signal 
(signal that is independent of isotropic dif¬ 
fusion, first term in the right hand of Eq. [1). 
The anisotropic signal within a voxel v is pre¬ 
dicted by demeaning the signal prediction in 
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Table 1: Mathematical notation and definitions for multiway arrays and decomposition. 


2 ^ A, w, b 


A tensor, a matrix, a vector and a scalar 


M-J.k), 

X(:,:,fc) 

X(„) G 


Entries of a tensor, a matrix and a vector 

Mode- 1 , mode -2 and mode -3 vectors are obtained by fixing all but one index 

Horizontal, lateral and frontal slices are obtained by fixing all but two indices 

Mode-n unfolding of multiway array X G whose entry at row i„ and 

column (ii - 1) J2 ■ ■ ■ In-iln+i ■ ■ ■ In H-h (/'n-i - l)fN + m is equal to 


Y = X x„ A G 


X G Xi Ai X2 A2 X3 A3 


X = i?ec(X) G 


So = diag{So{l),SQ{ 2 ),... ,Sq{K)) 


Multiway array by matrix product (in mode-Ji) where yi'i i'n “ 

Tucker decomposition: a 3 D multiway array X G is represented as the 

product of a core array G G by factor matrices A„ G 

Vectorization of multiway array X G with the entry at position 

h + EtL2[(4 - t)hh ■ ■ ■ 4 -i] equal to 

Diagonal Ny x Ny matrix with So(r^,r^) = So(r’) 


the following way: 

s{e,v)-h^J2zvfOf{e), (5) 

f€v 

where the mean is defined as 

k = ( 6 ) 

with Ng being the number of sensitization 
directions and Op(6) is orientation distribu¬ 
tion function specific to each fascicle, i.e. the 
anisotropic modulation of the diffusion sig¬ 
nal around its mean and it is defined as fol¬ 
lows: 

Of (8) = So{v) ^ . 

(7) 

The left-side hand in Eq. |5j the demeaned sig¬ 
nal, is the difference between the measured 
diffusion signal and the mean diffusion sig¬ 
nal. The right-hand side of Eq. is the pre¬ 
diction model. The LiEE model extends from 
the single voxel to all white-matter voxels in 
the following way: 

y w Mw, (8) 

where y G [g ^ vector containing the 

demeaned signal for all white-matter vox¬ 
els V and across all diffusion directions 9, 
i.e. ifi = S{6i,Vi) — Ivj- The matrix M G 
jfNgNyxNf coiuirui f the signal 

contribution Op(6) given by fascicle / at 
voxel V across all directions 6, i.e., M(i,/) = 


Op{6i), and w G contains the weights 
for each fascicle in the connectome. These 
weights are estimated by solving the follow¬ 
ing non-negative least-square constrained op¬ 
timization problem: 

min ^^||y “ Mw||^^ subject to wp > 0,V/. 

(9) 

This formulation of the LiEE model re¬ 
quires generating a matrix M (Eig. |2lE) 
that is very large and block-sparse. M is 
sparse because fascicles only cross a subset 
of all voxels. The matrix has NgNv rows by 
Np columns. Given an approximate human 
white-matter volume of 500 ml the number 
of voxels (Nv) can vary between 50,000 and 
500,000 depending on the spatial resolution 
of the dMRI acquisition (normally between 
4-1.25 mm^). Given that modern dMRI acqui¬ 
sition parameters measure between 30 and 
300 diffusion directions (Ng) and that whole- 
brain connectomes can contain as little as 
100,000 but as much as 10,000,000 fascicles 
(columns of M), the size of M can vary be¬ 
tween 10 and 100GB using standard double¬ 
precision floating-point format and sparse 
format. 

Eig. |2jD shows measurements of the size 
of M in GB for five individual brains and two 
datasets. Such large memory requirements 
necessitate large-memory compute-systems 
nowadays available in major research institu¬ 
tions (e.g., http://rt.uits.iu.edu/bigred2/, 
http://karst.uits.iu.edu). Below we intro- 
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duce a novel approach that dramatically re¬ 
duces the storage requirements of the LiFE 
model and makes LiFE suitable to run on 
standard desktop and notebook computers. 


11. Multiway decomposition of the 
Linear Fascicle Evaluation model 


Major contribution of the present work is 
a method to represent the LiEE model 
via Sparse Tucker Deco mposition (STD; 
iCaiafa and Cichocki 2012 )1. Before intro¬ 
ducing the method, we briefly discuss the 
standard and sparse Tucker-decomposition 
methods (see also Methods [T]|. 

Data compression by Tucker decomposition. 
Ledyard Tucker, a mathematician who spe¬ 
cialized in statistics and psychometrics, was 
first in proposing a decomposition approach 
to multidime nsional factorization problems 
i Tuckeilll966h . The approach provides a gen¬ 
eralization of the low-rank approximation 
for matrices to multiway arrays. Using the 
Tucker model, a 3-rd array X G IRkxLxL^ 
is approximated by the following decompo¬ 
sition: 


X w G Xi Ai X 2 A 2 X 3 A 3 , 


( 10 ) 


with a core array G G ]R,KixR2xR3 and factor 
matrices A„ G The decomposition is 

not guaranteed to provide a good aproxima- 
tion, but when it does it compresses data be¬ 
cause it uses a core array that is much smaller 
than the original multiway array, i.e. Rn <C 
I„. Indeed, the Tucker model provides us 
with a powerful compression method. This 
is because, instead of storing the whole orig¬ 
inal multiway array, we can st ore core array 
and factors. In this case, Eq. llOl i Tuck^ 1966 ) 
is called a low-rank Tucker model because G 
is small compared to X (see Eig. UK)- 

A simple example can help us explain 
how compression is obtained with a stan¬ 
dard low-rank Tucker model. Let consider a 
situation where R = = R 2 = R 3 and I = 

h = h = h- The low-rank Tucker decom¬ 
position of X requires storing only R^ -|- 3IR 
values. Where R^ is the number of entries of 
G and I and R the dimensions of Ai, A 2 and 
A 3 . Instead, the full multiway array would 
require storing values, compression is no¬ 


ticeable especially when G is small meaning 
that R 1. 



Figure 4: Classical and sparse Tucker dec omposi¬ 
tion. A. The classical Tucker decomposition iTuckeA 
It966l) allows representing a 3D multiway array X 6 
Rk X G X Ja product of core array (green) G 6 

RUXR 2 XR 3 \,y factor matrices A„ 6 RkxR„ 
low and blue). Data compression is achieved by consid¬ 
ering very small (dense) core arrays G, meaning that 
Rn <C L- B. The sparse T ucker Decomposition (STD; 
\Caiafa and Cichocki iZOlj) ). The core array G is large 
but sparse, not dense as in the classical Tucker model. 
Data compression is achieved because of the sparsity of 
the core array. 


Data compression by sparse Tucker decompo¬ 
sition: If the core G of a Tucker decomposi¬ 
tion of a multiway array is sparse, compres¬ 
sion c an be achieved even when G is very 
large jCaiafa and Cichockil. 12012 ). Consider 
a Tucker model with a large but sparse core 
array G G R^i XR2 xRs^ ggg p^g Only some 

entries in G are nonzero, i.e. 7 ^ 0 for 

m = 1,2,...M where M is the number of 
nonzero entries. Storing G all requires stor¬ 
ing: ( 1 ) the non-zero coefficients, ( 2 ) their 
location in Gj and (3) the factor matrices 
A„. This means that the storage cost (or¬ 
der) of the model is 4M -|- 3IR (assuming that 
R = R\ = R2 — R3 and I = f = I2 — L). 
Hereafter, we refer to this model as Sparse 
Tucker Decomposition (STD). Compared to 
the classical low-rank Tucker model, STD can 
provide better compression ratios with rela- 


8 










































C. F. Caiafa - F. Pestilli (May 2015) 


tively l ow M, i.e., with very spar se core ar¬ 
rays G dCaiafa and Cicho^. 2013h . 

Below we show how to apply the STD ap¬ 
proach to the LIFE model. To do so we first 
explain how the diffusion signals wifhin a 
voxel predicted by the LiFE model can be rep¬ 
resented using dictionaries of precompufed 
diffusion prediction signals. After that we 
show how to extend the decomposition ap¬ 
proach to the whole white matter volume. 


II.l Decomposing LiFE within a voxel us¬ 
ing prediction dictionaries 

The original LiFE model predicts anisotropic 
diffusion within a voxel, v, at each measured 
diffusion direction, 0, by using the orienta¬ 
tion of the white matter fascicles intersect¬ 


ing the voxel (Fig. |2p; IPestilli et all (l2014h l. 
Fascicles are defined as a lisf of {x,y,z) spa¬ 
tial coordinates within the brain, the fasci¬ 
cle nodes. To generate the prediction, a de¬ 
meaned tensor model is evaluated parallel 
to the orientation of each fascicle node (Eq. 
[7|. Hereaffer, we simplify fhe calculations by 
representing the signal prediction for any ar- 
bifrary fascicle-node orienfation using a dic¬ 
tionary of diffusion prediction signals gener- 
afed for a predefined sef of orientations regu¬ 
larly sampled over a grid. This grid covers a 
plausible rage of fascicles orientations on the 
unit-norm sphere (shown in Eig. |2]4) given 
the resolution of fhe dMRI dafa (Eig. 1^4). 

We define fhe dictionary matrix D G 
(Pig 1^) containing in its columns 
the demeaned canonical diffusion signals 
called atoms fhat can be used to approximate 
fascicles' contributions. More specifically, we 
define 


D(0,fl) = e 


_ ^-be‘Qae 


1 , 11 , 

Ne e 


where Qa is fhe diffusion prediction associ- 
afed fo the atom a, i.e. with a single ori¬ 
entation in 3D space. Dictionary atoms are 
identified by discretizing each spherical co¬ 
ordinate a and using a grid of values in 
fhe rage [0, n) (Eig. |^). The paramefer L in¬ 
dicafes fhe number of discretization samples 
on the grid. For example, L = 180 indicates 
a grid resolution of 1° x 1° wifh an approxi- 


mafe atom number Nn ~ L^. 

The dictionary D G Ji^JgxNa allows pre¬ 
dicting the matrix Mj, G defined as 

a block of matrix M corresponding to voxel 
V (see Fig. |5j2), by decomposing it in the fol¬ 
lowing way: 


= So{v)-D^^, (12) 

where G is a sparse matrix whose 

non-zero entries at column / indicates the 
dictionary atoms selected to predict the voxel 
signal, given the orientations of fhe fascicles 
in fhe voxel. In sum, instead of sforing the in¬ 
dividual contribution of each fascicle wifhin 
a voxel, we store only the indices to the dic¬ 
tionary atoms, and avoid multiple versions 
of very similar signals. Below we exfend the 
STD of LiFE from single voxels to the entire 
white matter volume. We show that for each 
dafa sef used we can find a finite number of 
dictionary atoms that generate predictions as 
accurate as those obtained with the original 
LiEE model. 


11.2 Extending the LiFE decomposition 
model across all voxels 


The original LiFE model comprises a single 
very large block-sparse matrix (M in Eq. |^. 
The matrix M G can be converted 

into a multiway array by mapping the dif¬ 
fusion directions (0), voxels (v) and fasci¬ 
cles (/), onto the dimensions of a 3D mul¬ 


tiway array M G 


X Nt, X N, 




Representing 


M as a 3D multiway array M allows us us¬ 
ing a Sparse Tucker dec omposition approach 
and compress its size dCaiafa and Cichocki , 

201 2h . 


To fully exploit the Sparse Tucker Decom¬ 
position we consider both the LiEE model M 
and its optimization problem, as defined in 
Eq. 1^ This equation can be conveniently 
rewritten using multiway arrays in the fol¬ 
lowing marmer (see Eig. [^): 


Yf«Mx3 W^ (13) 


where matrix Y G is the matrix ver¬ 
sion of vector y, and M G jg ^ g. 

way array. The laferal slices of M are defined 
by fhe matrices G defined in Eq. 

[121 these represent matrices correspond- 
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A 



B 



Ng 


II 



Demeaned predicted diffusion signal 
for three example atoms 



IS 


Nf 

Diffusion signals associated to fascicles 
intersecting voxel V 


□ Empty matrix entries 



□ □ Nodes of fascicle/i within voxel ft 

□ Node of fascicle /2 within voxel V 

□ Empty matrix entries 


Figure 5: Discretization of the LiFE model. A. Discretization of the unit-norm sphere and mapping of the spher¬ 
ical coordinates to a single index a. The unit sphere of diffusion directions is sampled by using a uniform grid 
of spherical coordinates (using a selected choice of a and f values see Fig. |2]4 for definitions). B. Construction 
of the diffusion prediction dictionary. Demeaned diffusion predictions are generated for each spherical coordinate 
(orientation) and measured diffusion direction. Predictions are stored in the columns of a dictionary of diffusion 
predictions, D, whose columns identify different fascicles orientations (Na) and rows the predicted diffusion in a 
measured diffusion direction (Ng). C. Modeling the diffusion signal in a voxel using the discretized LiFE model. 
The measured signal in voxel v is organized as a matrix, Mv. The signal predicted by each fascicle (two in our 
example) and nodes (three in our example) passing in v is approximated by combining the diffusion prediction of 
the dictionary atoms (columns ofD) with directions closest to the orientation of the fascicles nodes (yellow, blue 
and red). Non-zero entries in 4>t, indicate the atoms corresponding to the nodes in the fascicles, fi and f 2 in the 
example. 


ing to each white matter voxel. Furthermore, 
Eq. [12] shows that each lateral slice of M is 
scaled by So(i^), which we represent at a di¬ 
agonal matrix Sq. This allows us writing the 
LiFE model that encompasses all voxels by 
using an efficient Sparse Tucker Decomposi¬ 
tion (STD): 

M = ^XiDX2 So, (14) 

where the 3 -way array O G -^NaxNvxNf 
has as lateral slices the matrices 
Oj, G matrix Sq = 

diag{So{l),So{2),...,So{N,)) G is a 

diagonal matrix with values So(u) along the 
main diagonal. By combining Eqs. (IT3ll and 
ifTH we obtain the following Sparse Tucker 
decomposition of the LiEE model (see also 


Fig- 1^): 

Y « O Xi D X 2 So X 3 w^. (15) 

is the core of the multiway decomposition, 
because O is sparse it results in strong data 
compression (see section |Vj. Below we intro¬ 
duce the set of operations necessary to build 
and optimize (fit) the STD LiFE model. 

III. Building the STD LiFE model 

The STD LiEE model is built using a large 
collection of white matter fascicles estimated 
using co mputational t r actog raphy, the con- 
nectome dSporns et a i |200^ . Eascicles in 
the connectome are represented as a list of 
(x, y, z) brain coordinates. The O core-array 
of STD LiEE model is then built in the follow¬ 
ing way: ( 1 ) the orientation of each fascicle's 
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(/) node (n) is identified in spherical coordi¬ 
nates (a, /3); (2) nodes' orientation is mapped 
to the closest atom (a) in the dictionary (D); 
(3) the entry in O corresponding to each iden¬ 
tified voxel, fascicle and atom is set to 1: 

0(fl,i;,/) = l, (16) 

The rest of the entries in 4> are set to zero. 



Figure 6: Sparse Tucker Decomposition of LiFE prob¬ 
lem. A. Multiway version of the LiFE method y k, 
Mw, Eq. [I4] B. LiFE problem decomposition by 
Sparse Tucker method, Eq. HSl 


where M e u^e^vxNf original LiFE 

model, w e the fascicle weights and 
y G R^e^D the demeaned diffusion signal. 
Because the STD version does not explic¬ 
itly store M in the following section we de¬ 
scribe how to perform two basic operations 
(y = Mw and w = M^y) using the multiway 
decomposition to compute the optimization 
gradient. Appendix A reports pseudocode 
implementing the operations. 


IV.l Computing y = Mw 

The product Mw can be computed in the fol¬ 
lowing way using a 3D-array by vector prod¬ 
uct: 

Y = Mx 3 W^ (18) 

where the result is a matrix Y G ^ 

matrix version of the vector y. Using the STD 
LiFE model (Eq. [T4ll the product is written as 
follows: 

Y = ^ Xi D X 2 So X 3 w^. (19) 


IV. Optimizing the STD LiFE model Computing w My 


Building the STD LiFE model is the first 
of two steps in evaluating a brain con- 
nectome. The final step requires finding 
fhe non-negative weights that least-square 
fit the measured diffusion data. This is 
a convex problem that can be solved us¬ 
ing a variety of Non-Negative Least Squares 
(NNLS) optimization algorithms Q. The 
original LiFE problem, was solved using a 
NNL S algorith r n bas ed on first-order meth¬ 
ods (I Kim et all, l2013h . Hereafter, we show 
how to modify fhe opfimizafion algorifhm to 
exploit the STD LiFE model. 

The optimization of the STD LiFE needs 
to be performed using its core array (O) and 
matrices (D and Sq; Eq. iMll . The gradient 
of the original objective function for the LiFE 
model can be written as follows: 

Vw Q||y - Mw|M = M^Mw, (17) 


The product w = M^y can be computed us- 
ing the STD LiFE model in the following way 
( Kolda and Bader , 2009 ): 


w = M^y = M( 3 )y = (Sq ® D^)y, ( 20 ) 

where 0 is the Kronecker product. Eq. 
equation can be written as follows: 

w = «I>( 3 )Z;ec(D’’YSo), ( 21 ) 


where vec{) stands for the vectorization op¬ 
eration, i.e. to convert a matrix to a vec- 
tor by stacking its column s in a long vector 
( Caiafa and Cichocki . 2012 ). 


Because matrix is very sparse, we 
avoid computing the large and dense matrix 
D^YSq and multiply only the non-zero en¬ 
tries in ‘1*(3). This allows maintaining effi¬ 
cient memory usage and limits the necessary 
number of CPU cycles. 
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Figure 7: Reduction in memory usage using STD 
LiFE. A. The memory requirement o/M plotted against 
memory requirement of the STD model (Eq. [I^ L = 
360, STN96 data, probabilistic tractography, Lmax = 
10). B. Memory requirement for M and STD (Eq. [I^ 
as function of the number of measured diffusion direc¬ 
tions (FLiheta^ T = 360, STN96 data, probabilistic trac¬ 
tography, Lmax = 10). C. Memory requirement for M 
and STD (Eq. \1^ as function of the number of fas¬ 
cicles in the connectome (Nf, L = 360, STN96 data, 
probabilistic tractography, Lmax = 10) 


M and the STD model can also be computed 
analytically. To do so we assume that all fas¬ 
cicles have the same number of nodes Nn and 
that there are no more than one node per fas¬ 
cicle, per voxel. Under these assumptions the 
amount of memory necessary to store each 
fascicle / is proportional to NgNn, thus the 
storage cost of M is: 

C(M) = 0{NnNgNf). (22) 


Conversely, storing fascicles in the STD 
model require 4Nn values only (i.e. the set of 
the non-zero coefficients and their locations 
within the core multi-way array O). Thus 
the amount of memory required by the STD 
model is: 


C(M) = OiiNnNf + NgNa), (23) 


where NgNa is the storage associated with 
the dictionary matrix D € please 

note that Sq is absorbed as O = ® X 2 Sq 
without affecting sparsity, i.e., with no ef¬ 
fect on the storage. Storage reduction can be 
straightforwardly computed as follows: 



Ng 

N„Nf- 


(24) 


Fig . [ZB shows rnemor y usage for the mod¬ 
els in ( Pestilli et al , 2014 ) and Eq. [14] as func¬ 
tion of number of diffusion directions Ng in 
the data. Given a fixed number of fascicles 
Nf and nodes Nn in a brain storage of M 
grows linearly with the number of diffusion 
directions see Eq. |22j much faster 

than the STD model (Eq. l23l . 

Eig. 13 ^ shows memory usage for the 
models in llPestilli et al , 2014 ) and Eq. |14| as 
function of the number of fascicles in a con¬ 
nectome (Nf). The storage of M grows lin¬ 
early with Nf and grows much faster than 
the STD model. 


V. Model storage reduction 


Fig. 13^ compares mem ory usage by M (Fig. 
Et; IPestilli et all (l20l4 ) and the STD LiFE 
model (Eq. (Till . The storage requirements of 


In sum, the STD model provides substan¬ 
tial reduction in memory requirements. The 
reduction in memory consumption achieved 
by the STD model becomes more important 
as the number of measured directions and 
fascicles in a connecfome increase. 
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VI. Model accuracy 

The STD model provides only an approxima¬ 
tion of the original LIFE model. This is be¬ 
cause of the discretization introduced by the 
dictionary (D, see Fig. |5j\). Fig. shows that 
the STD model achieves accuracy similar to 
the original model. We compared the accu¬ 
racy of the STD model in both, predicting the 
demeaned signal and estimating the weights 
of the original model. 



Figure 8: The STD model is as accurate as the origi¬ 
nal LIFE model. A, Difference in mean cross-validated 
r.m.s error for a range of discretization steps between 
the STD and original LiFE models. B. Relative error in 
estimating the fascicles weights between the STD and 
original LiFE model for a range of discretization steps 
(STN96 data, probabilistic tractography, Lmax = 101. 


Fig. |8j4 shows the difference between 
the r.m.s error of the STD and original LiFE 
model. We fit the STD model using dif¬ 
ferent number of discretization steps of the 
spherical coordinates (l5]4). To do so, the 
number of discretization steps was varied 


(L = 23,45,90,180,360 and 720) and the dif¬ 
ference in mean r.m.s. in predicting the de¬ 
meaned diffusion signal across the whole 
white-matter volume was computed in refer¬ 
ence to the mean r.m.s. of the original model. 
The mean r.m.s. used for comparison was 
c ross-validated to an independent data set 
( dPestilli et ail2014h l. 

Eig. 1^ shows the relative error of the 
STD model in estimating the weights as¬ 
signed by the original LiEE model to each 
fascicle. We fit the STD model using dif¬ 
ferent number of discretization steps of the 
spherical coordinates as in|8j4. We computed 
the relative error in the estimated weights 
as IIwg — wi||/||wo|| where wg, are the 
weights of the original and STD models, re¬ 
spectively. Results show that for the data sets 
tested, L > 360 allows achieving error below 
1 %. 


VII. Reproductio n of results from 


Pestilli et all (120141) 


We demonstrated that a Sparse Tucker De¬ 
composition model achieves similar accuracy 
to the original LiEE model and reduces stor¬ 
age requirements by 97%. Eig. shows 
that the STD mod el replicates a ma jor re¬ 
sult of the original llPestilli et al l2014h . The 
scatter plot shows the cross-validated r.m.s. 
in predicting the demeaned diffusion signal 
of a probabilistic and deterministic tractogra¬ 
phy cormectome. Results show t hat the STD 
mode l replicates the findings of Pestilli et all 
(12014 ) demonstrating a larger r.m.s. for the 
deterministic tractogrpahy cormectome in a 
majority of the white matter volume (see Fig. 

for the same plot computed using the 
original LiFE model). 


VIII. Open source LiFE software 
and reproducibility of results. 

The Matlab implementation of LiEE soft¬ 
ware using the new STD model is pro¬ 
vided at github.com/brain-life/life and 
francopestilli .github. io/life The Matlab 
impl ementation uses th e Matlab Tensor Tool¬ 
box ( Bader et~^. 12012 ) and mex files com¬ 
piled for Mac OSX and 64-bit Linux dis¬ 
tributions. The software has been tested 
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to run on a standard Notebook computers 
with less then 8GB of RAM. Processing a 
whole brain connectome with 500,000 fasci¬ 
cles of the STN150 dataset using a single 
CPU thread (Notebook computer, 2.2GHz In¬ 
tel Core i7 processor and 8GB RAM) requires 
about 5 hours (L = 360). 



10 50 90 

Probabilistic 


Figure 9: The STD model replicates major results of 
jPestilli et aj \20l4) . The scatter plot shows the r.m.s. 
errors in predicting the demeaned diffusion signal of 
two tractography models. The r.m.s. error of a proba¬ 
bilistic connectome is plotted against that of a determin¬ 
istic connectome. R.m.s. was computed using a single 
subjects from the STN96 data set using CSD-based 
probabilistic tractography (Lmax = 101 and tensor- 
based deterministic tractography and the STD LiFE 
model with discretization parameter L = 360. 


IV. Discussion 


The number of brain dataset collected using 
modern neuroimaging methods is growing 
at exponentially fast pace. At the same time 
both data resolution, as well as the size of the 
populations of human br ains being acquired 
and share d are gr ow ing dAmunts et ai . 2013 
Scott et al 2011: _ Van Essen and Ueurbil 


2012 : Van Essen et al . 20131: Zuo et ah 2014 ). 

The next generation of brain science will 
require collaborative efforts between the neu¬ 
roscience community as well as the infor- 
matics and signa l proce s sing communities 
I Garvfallidis et all. 12014 iGoreolewski et al , 
2 OIII: Perez and Grangei , 200?!) . 

The present article focuses on the appli¬ 
cation of modern signal processing meth¬ 
ods for big data to neuroimaging. More 
speficically we present an example appli¬ 
cation to diffusion imaging and tractogra¬ 


phy evaluation ) Pestilli et all, 12014 ). Ma¬ 


jor value of dMRl and tractography is 
to allow measuring white-matter in living 
brains in an individualized manner, one 
brain at the time. Major efforts are be¬ 
ing put forth to improve the representation 


and micro-structural level dAssaf and Bassei 

2OO.5I: 1 Assaf et al. 

2008, 

201 .31: 1 

laducci et al 

2015allbl: iGarvfallidis et a 

II2OI2I: 

Pestilli et al 

20141 Smith et al 1 

2006: 

Yeatman et all 

2012 

Yendiki et al 2011 

Zhang et al 

2012) 

The 


digital nature of neuroimaging data and the 
availability of large in-vivo databases affords 
developing new methods to built individ¬ 
ual connectomes and co mpute their accura cy, 
and validate the results jPestilli et al . 2014 ). 


Computational tractography methods ex¬ 
ploit the diffusion-weighted signal measured 
with MRl to identify plausible trajectories 
of white matter tracts. Most tractogra¬ 
phy methods estimate candidate tracts one 
node at a time. Beyond the great excite¬ 
ment brought about by all these technolo¬ 
gies, much work is necessary to develop ap- 
proache s for evaluation and validation of 
results d Catani et^. 2012L lones et 12013 : 
Sporns et al 2005 ). To date two primary ap¬ 


proaches have been used to tractography val¬ 
idation. First, the geometric accuracy of trac¬ 
tography has been comp ared to simulated 
and physical phantoms llFillard et al , 2011 


ISchreiber et al , 20l4lzheng et ai 2014 ). Sec 
ond, tracts and connections identified us¬ 
ing tractography have been compared to 
connections obtained with staining meth¬ 
ods i n postmortem tissue jAzadbakht et al , 
201 . 5 !: ^rker et ai 20021: Seehaus et al . 201.3 


Sherbondy et ai 120081: Thomas et all. 20141) . 
These validation approaches have helped es¬ 
tablishing that tractography is accurate to a 
certain degree and can well identify tracts 
within the core white matter. 

Attempts to improve tractography have 
focussed on what is called global tractog¬ 
raphy. Global tractography methods eval¬ 
uate the plausibility of the tractography 
estimates by comparing individual stream¬ 
lines trajectories with models of the dMRl 
signal or heuristic rul es abo u t streamline 
smoothness dAganj et al , 201 ll: Fillard et ai 
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2009, 

201 it Ibabdietal, 

20071 20081; Lietal, 

2012 : 

Neher et al. 

2012 : 

Reiser! et al. 

2011 

Sherbondv et al. 

200 sl). 

Alternative ap- 


proaches to global tractography have pro¬ 
posed the concurrent estimation of the bio¬ 
physical properties of the whit e matter fas¬ 
cicles and surroun d ing ti s sue dCirard et al , 


2015 

: Kadenetal, 

20151: 

Sherbondv et al 

2011 

, 20091: Izhang 

and Laidlawl 

2006). His- 


torically global tractography methods have 
been considered computationally intensive 
tasks with limited applicability to routine 
study of fhe human brain in healfhy and dis¬ 
eased populafions. 


Recently a method for tractography 
evaluation b ased on linearized models 
llPestilli et all 2014 ) has been proposed 
and applied t o stu d y living conne ct omes 


I Gomez et all, [201^ Takemura et al , 20151: 
Yeatman et at l2014h . The method separates 
the process of fracking fro m that of global 
evalu ation of fhe fascicles dSherbondy ef al 
2008 ). Similarly a family of algorithms 
based on linearized models have been pro¬ 
posed to solve a series of estimation prob¬ 
lems identified in modern fractogrpahy, as 
well as for studying of fhe tissue microstruc¬ 
ture ( Daducci et al l2015albl) . This new gen¬ 
eration of linearized methods for microstruc¬ 
ture estimation, tractography evaluation are 
paving the road for new avenues of rnvesti- 
gafion and fhe study of fhe whife mafter in 


VIVO. 


The contribution of the work presented 
here is to introduce the general framework 
of multiway decomposition methods that 
can potentially be applied to any type of 
linearized neuroimaging models. The de¬ 
composition methods allow reducing the 
computational complexity of the neuroimag¬ 
ing models and in doing so can increase the 
impact of these models. The application of 
multiway decomposition methods to neu¬ 
roimaging will pave the road to use global 
estimation and evaluation methods to study 
large populations of human brains wifh in¬ 
creasingly high spatial resolution and to ap¬ 
ply the modern methods to study the hu¬ 
man brain in normal and clinical populations 
such as tho se avail a ble in modern databases 
iScott et al , 2011 : Van Essen and Ugurbill. 


2OI2I: Van Essen et al , 2013 ). 


Appendix A: Computational 

ALGORITHMS. 

Below we report pseudo code for fhe fwo 
operations necessary to fit the decomposed 
LiFE model. These operations provide the 
option of being implemenfed using multi¬ 
threading methods. This is because the in¬ 
stances in the loop are independent and al¬ 
low parallel computations. 


Algorithm 1: y = M_times_w(^,D,So,w) 

Require: Decomposition components D, So) and 
vector w 6 , 

Ensure: y = Mw 

1: O = $ X 2 So; the result is a very large but 
still very sparse 3D-array. 

2: Y = $ X 3 w^; the result is a large but very 
sparse matrix (N^ X Nb) 

3: Y = DY; the result is a relatively small 
matrix (Ng X Ny) 

4: y = hec(Y) 

5: return y; 


Algorithm 2 : w = Mtransp times ylO.D.Sn.y) 

Require: Decomposition components D, So) and 
vector y e 

Ensure: w = M^y 

1: Y 6 t— y 6 R^o^”; reshape vector y 

into a matrix Y 

2: O = ^ X 2 So ; the result is a very large but 
still very sparse 3D-array. 

3: [a, V, f, c] = get_nonzero_entries(^); a(n),v(n), f(n), 
c{n) indicate the atom, the voxel, the fascicle and 
coefficient associated to node n, respectively, with 
n = 1,2,...,N„; 

4: w = 0 6 R^f; Initialize weights with zeros 

5: for h = 1 to N„ do 

6: w{f{n)) = w(f{n)) + D'^(:,a(n))y(:,v(n))c(n); 

7: end for 

8: return w; 
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